pacman::p_load(sf, spNetwork, tmap, tidyverse)Advanced SPPA Methods - NCKDE
Getting Started
Data Sources
Data for this exercise are from public sources and will be used to analyse the distribution of childcare centers in the Punggol planning area. Two datasets in ESRI shapefile format will be used:
A line feature geospatial dataset which includes the road network of Punggol planning area
A point feature geospatial dataset which includes the location of childcare centers in the Punggol planning area
Installing and launching R packages
This exercise will make use of five R packages: sf, spNetwork, tidyverse, and tmap.
sf - for importing, managing and processing vector-based geospatial data
tidyverse - collection of packages for performing data importation, wrangling and visualization
tmap - for plotting cartographic quality maps
sPNetwork - provides functions for performing SPPA methods like KDE and K-function on a network. The package can also be used to build spatial matrices to conduct traditional spatial analyses with spatial weights based on reticular distances
The code chunk below uses p_load() of pacman package to check if the packages are installed in the computer. It installs them first if they are not. It then loads them into R.
We use the random seed 1234 to ensure reproducibility of results
set.seed(1234)Data Import and Preparation
The code chunks below uses st_read() of the sf package to load the street and childcare data into their respective dataframes.
network <- st_read(dsn="data/geospatial",
layer="Punggol_St")Reading layer `Punggol_St' from data source
`C:\drkrodriguez\ISSS626-GAA\In-class\In-class_Ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TMchildcare <- st_read(dsn="data/geospatial",
layer="Punggol_CC")Reading layer `Punggol_CC' from data source
`C:\drkrodriguez\ISSS626-GAA\In-class\In-class_Ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TMWe examine the contents of these objects by running the code chunks below
networkSimple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
First 10 features:
LINK_ID ST_NAME geometry
1 116130894 PUNGGOL RD LINESTRING (36546.89 44574....
2 116130897 PONGGOL TWENTY-FOURTH AVE LINESTRING (36546.89 44574....
3 116130901 PONGGOL SEVENTEENTH AVE LINESTRING (36012.73 44154....
4 116130902 PONGGOL SEVENTEENTH AVE LINESTRING (36062.81 44197....
5 116130907 PUNGGOL CENTRAL LINESTRING (36131.85 42755....
6 116130908 PUNGGOL RD LINESTRING (36112.93 42752....
7 116130909 PUNGGOL CENTRAL LINESTRING (36127.4 42744.5...
8 116130910 PUNGGOL FLD LINESTRING (35994.98 42428....
9 116130911 PUNGGOL FLD LINESTRING (35984.97 42407....
10 116130912 EDGEFIELD PLNS LINESTRING (36200.87 42219....childcareSimple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
First 10 features:
Name geometry
1 kml_10 POINT Z (36173.81 42550.33 0)
2 kml_99 POINT Z (36479.56 42405.21 0)
3 kml_100 POINT Z (36618.72 41989.13 0)
4 kml_101 POINT Z (36285.37 42261.42 0)
5 kml_122 POINT Z (35414.54 42625.1 0)
6 kml_161 POINT Z (36545.16 42580.09 0)
7 kml_172 POINT Z (35289.44 44083.57 0)
8 kml_188 POINT Z (36520.56 42844.74 0)
9 kml_205 POINT Z (36924.01 41503.6 0)
10 kml_222 POINT Z (37141.76 42326.36 0)We see from the output of calling childcare that the dataset includes a third dimension– a z-coordinate. The spNetwork package can only accept objects in (x,y) and not in (x,y,z)
To go around this, we re-import the data with st_zm() (using default arguments: drop = TRUE, what = “ZM”) to remove the z-coordinate
childcare <- st_read(dsn="data/geospatial",
layer="Punggol_CC") %>%
st_zm()Reading layer `Punggol_CC' from data source
`C:\drkrodriguez\ISSS626-GAA\In-class\In-class_Ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TMChecking the content of the object shows that the third dimension is now dropped
childcareSimple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
Projected CRS: SVY21 / Singapore TM
First 10 features:
Name geometry
1 kml_10 POINT (36173.81 42550.33)
2 kml_99 POINT (36479.56 42405.21)
3 kml_100 POINT (36618.72 41989.13)
4 kml_101 POINT (36285.37 42261.42)
5 kml_122 POINT (35414.54 42625.1)
6 kml_161 POINT (36545.16 42580.09)
7 kml_172 POINT (35289.44 44083.57)
8 kml_188 POINT (36520.56 42844.74)
9 kml_205 POINT (36924.01 41503.6)
10 kml_222 POINT (37141.76 42326.36)Visualization of the sf objects
The code below uses plot() in one map. The add=T argument in the second line allows the two plots to be added one over the other.
plot(st_geometry(network))
plot(childcare, add=T, col="red", pch=10)
The following code chunk produces a similar map using tmap package.
tmap_mode('view')tmap mode set to interactive viewingtm_shape(network) +
tm_lines() +
tm_shape(childcare) +
tm_dots(col = "red")tmap_mode('plot')tmap mode set to plottingNetwork Constrained Spatial Point Analysis
Preparing the lixels objects
A requirement for NKDE is that the lines objects needs to be cut into lixels. The code below uses lixelize_lines() on network, using a length of 700 and a minimum distance (mindist) of 350.
lixels <- lixelize_lines(network,
700,
mindist = 350)Generating the line center points
The code below generates the center points for the lixels
samples <- lines_center(lixels)xxx
tmap_mode("view")tmap mode set to interactive viewingtm_shape(lixels) +
tm_lines() +
tm_shape(samples) +
tm_dots(size = 0.01)tmap_mode("plot")tmap mode set to plottingComputing NKDE
The code below computes for the NKDE of childcare centres around network
densities <- nkde(network,
events = childcare,
w = rep(1, nrow(childcare)),
samples = samples,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5,
sparse = TRUE,
verbose = FALSE)The output of the code is a list of density or intensity values. We copy these values into the original lixel and lixel midpoint dataframes.
samples$density <- densities
lixels$density <- densitiesNote that the previous values are based on metres and resulted in very low density values. We can change the density values to event per square km by using the code below
samples$density <- samples$density*1000
lixels$density <- lixels$density*1000Computing the K and G Functions
The code block below computes for the K- and G-functions based on the data using kfunctions()
kfun_childcare <- kfunctions(network,
childcare,
start = 0,
end = 1000,
step = 50,
width = 50,
nsim = 49,
resolution = 50,
verbose = FALSE,
conf_int = 0.05)The output can be visualized by calling plotk for the K-function and plotg for the G-function
kfun_childcare$plotk